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ABSTRACT 


The search for more efficient gas turbine engines 
has led to the proposal for the replacement of metal 
high temperature components with ceramic components. 
Essential to this effort is the numerical analysis of 
proposed designs. This thesis report describes the model 
discretization of a proposed blade design, the development 
of pre- and post-processors for the ADINA finite element 


code and the initial stress analysis of the modeled blade. 


OISTRIBUTION /AVAILABILITY GODES 


_Uitt.__ AVAIL ender SPECIAL 


a 


' 
a oe ee ae 


1 BES 


IV. 


4 4 


TABLE OF CONTENTS 


INTRODUCTION ---~-------------------------~----- 
ANALYTICAL MODEL DEVELOPMENT ------~------------ 
A. BLADE NOMENCLATURE ~-----------~------------- 
B. GEOMETRIC DEFINITION -----------------~----- 
1. Airfoil ~-------~-------------~----~----- 
2. Fillet --------~-------+-----------~----- 


a. Determination of the External 


Normal ----~-- 


b. Locus of 0.3 Inch Radius Centers --- 


ec. Locus of 0.1 Inch Radius Centers --- 


d. Definition of Perimeter Point at 


Arbitrary Z-Level ~------~---------- 
3. Attachment Root ------~----------------- 
C. ANALYSIS MESH DISCRETIZATION -----~--------- 


PRE-PROCESSOR DEVELOPMENT 
A. GENERAL DESCRIPTION OF 
B. MODIFICATIONS TO ADINA 
C. CENTRIFUGAL LOADING -- 

1. Mathematical Formu 


2. Computer Algorithm 


Centrifugal Load - 

a. Input ~-<-=--- 

b. Calculation of 

Loads +------~ 

c. Output ------- 

PROBLEM SOLUTION --------- 


POST=-PROCESSOR DEVELOPMENT 


ADINA ----------..--- 
tation ==—s=<«sssceace 


for Program 


Ce ieee ee ee et 


10 
13 
13 
14 
14 


14 


15 
15 


16 


16 
16 
alg) 
27 
ZT 
29 
30 


32 


33 


34 


34 
34 
37 
38 


ahh a Ne 


SS ON Re 


PROGRAM MOCOW? cneeedeeneeneserecensseacecees 


PROGRAM GERESG o2acanveensencecesacaauese ous 


PROGRAM CONTOUR PLOT DATA ------------------ 


PROGRAM CONTOUR PLOT <2<oco~aeeeeecnconcance 


VI. RESULTS OF ADINA STRESS ANALYSIS --------------- 


A 
B. 


COMPARISON OF MAXIMUM PRINCIPAL STRESSES --- 


CONTOUR PLOT ANALYSTS =--<--2----.---.-5-5. 


WET) CONCLUSIONS =2=2262 39a een oe eee ee 


A. 
B. 


APPENDIX 


APPENDIX 


APPENDIX 
APPENDIX 
APPENDIX 


APPENDIX 


AMAINSTS RESULTS. 22os aa i ee eee 


OPPORTUNITIES FOR FURTHER RESEARCH --------- 


A 


Cc 
D 
E 


E 


ADDENDUM TO ADINA USER'S MANUAL, REPORT 
82448-1, MIT, SEPTEMBER 1975 ------------ 


CONVERSION OF CALCOMP PLOTTING ROUTINES 
FOR USE ON A TEKTRONIX 4012 TERMINAL ---- 


PROGRAM CENTRIFUGAL LOAD LISTING -------- 
PROGRAM KCONT LISTING --<--<.-..-.~<- <= 
PROGRAM STRESS LISTING «oo-+e-2--4--seeuu 


PROGRAM CONTOUR PLOT DATA LISTING ------- 


LIST OF REFERENCES -2ccksuesesew ween scemwesecueunnwse 


INITIAL DISTRIBUTION LIST ooe<encencwceseenceecn snes 


39 


40 


LIST OF TABLES 


1. ADINA PROBLEM SOLUTION TIMES 


2. COMPARISON OF RESULTS OF ANALYSES USING TWO, 
THREE AND FOUR POINT INTEGRATION 


LO. 


ll. 


LIST OF FIGURES 


Illustration of Nomenclature and Reference System 
of Turbine Rotor Blade -------~----------~------- 


Illustration of Plane in Which Fillet Geometry 
is Defined 


Page 1 of 2 ---------~---------- = aera ee ee 
Page 2 of 2 ---------~----------~---------------- 
Illustration of Plane in Which Attachment Root 

Geometry is Defined -~--------------------------- 
Airfoil Mesh, Elements 1 through 12 -------------- 


Attachment Root Mesh 


Page 1 of 2 ~-----------+----------~--------~+----- 


Page 2 Of 2 ~--------- 9-9-2252 93-95 ------------ 


Completely Assembled Finite Element Analysis 
Mesh of Proposed Ceramic Gas Turbine Blade 


Design -------------- 92-52-25 3 -- = - = $+ -- 


Stress Output Locations for ADINA Linear Static 


Analysis Results ------~-----------~-------------- 


Region of Highest Stress Concentration, Elements 


60, 61 and 62 ~--------~--+--------~----------~--- 


Attachment Root Finite Element Mesh ---------~+--- 


Mid-Plane Dovetail Cross-Section 


Page 1 of 3 -------------~+----------~------------ 
Page 2 of 3 -----+--------~--+--------~------------ 


Page 3 of 3 ----~--------~---------- += -- = -- =~ ++ 


Dovetail Pressure Side 


Page 1 Of 3 ----- 2020-4202 -5-5 5-05-25 ------- = +--+ 
Page 2 Of 3 --------- 22-23-22 5 oo ee ++ 


Page 3 Of 3 a<---eensennewen nee en en wn ewe sene---=- 


19 


20 


21 


22 


23 


24 


25 


26 


36 


47 


48 


49 
50 


$1 


52 
93 
54 


12. Dovetail Suction Side 


Page 1 of 3 ----------~------------------------- 55 
Page 2 of 3 ----------9--- 2-2-2 22 - eo = 56 
Page 3 of 3 ----------~------------------------- 57 { 


13. Level 2.8404 


Page 1 of 3 ----------~------------------------- 58 
Page 2 of 3 ----------+---------4--------------- 59 
Page 3 of 3 ----------~------------------------- 60 


14. Level 2.9029 


Page 1 of 3 ------------------------------------ 61 
Page 2 of 3 -----------------2------------------ 62 
Page 3 of 23 ------------------------------------ 63 


15. Level 2.9654 


Page 2 of 3 -----------------------+------------- 65 
Page 3 Of 3 ------------- 22-2222 += == === 66 


Page 1 of 3 NE fee eae 22-22 -- +--+ ------ +--+ 64 
| 
| 


16. Maximum Principal Stress Contours for Pressure 
Side of Airfoil (ksi) -------------------------- 67 


17. Maximum Principal Stress Contours for Suction 
Side of Airfoil (ksi) -------------~------------- 68 


i 


ACKNOWLEDGEMENTS 


The author wishes to express his appreciation to the 
faculty and staff of the Mechanical Engineering Department 
of the Naval Postgraduate School for their support and 
motivation throughout his work at this institution. In 


particular, a great debt of gratitude is owed to Dr. Gilles 


Cantin, Professor of Mechanical Engineering, for his guidance 


and friendship as instructor and thesis advisor. Thanks are 
also extended to the staff of the W. R. Church Computer 
Center of the Naval Postgraduate School for their invaluable 
assistance. 

Finally, the author thanks his wife, Shay, and daughter, 
Heather, for their support and forebearance during the 


demanding period of his academic pursuit. 


I. INTRODUCTION 


The need to increase fuel efficiency of gas turbines has 
led to a sizable research investment into suitable materials 
and designs of ceramic gas turbine components to replace 
high temperature resistant superalloys. This replacement is 
envisioned to yield a threefold benefit: 

1) The refractory characteristics of ceramics will allow 
a greater source-sink temperature differential increasing 
plant efficiency without penalizing overall performance by 
requiring increased cooling air bled off the compressor. 

2) The lower specific weight of ceramics relative to 
the metal components to be replaced will lower plant weight 
thereby increasing the plant's power-to-weight ratio. 

3) Ceramic constituent components are, in general, 
domestically available giving them a strategic edge over 
superalloys. 

The basic problem preventing realization of these benefits 
is the brittle characteristic of ceramics. Since the brittle 
property is associated with the desirable refractory char- 
acteristic, material research and development can only be 
expected to yield marginally better ceramics in the near 
future. Any short term success will, therefore, be a result 
of design innovations which will compensate for the brittle 
nature of ceramics to allow their use for components which 
have formerly been engineered on the basis of use of ductile 


materials. 


10 


Essential to this design effort is the use of numerical 
techniques to model proposed designs and develop alternative 
criteria with the ultimate goal of minimizing stress concen- 
trations which key catastrophic brittle failure. This thesis 
describes the efforts of the author to utilize the computer 
hardware and software available at the Naval Postgraduate 
School (NPS) to analyze stress distributions in gas turbine 
ceramic replacement components. Time constraints dictated 
limiting the scope of analysis to the modeling and develop- 
ment of pre- and post-processors for a proposed first stage 
ceramic turbine blade subjected to centrifugal loading. 

A turbine blade design, developed by the AIRESEARCH 
division of GARRETT CORPORATION under a Navy managed contract 
for development of a one-hundred hour test gas turbine 
engine using ceramic components for the combustor, stator 
vanes, and rotor blades, was used in this work. Stress 
distribution analyses were centered around the "Automatic 
Dynamic Incremental Nonlinear Analysis" (ADINA) code devel- 
oped by Dr. Klaus-Jurgen Bathe of the Massachusetts Institute 
of Technology [Ref. 1] as implemented on the NPS IBM 360-67 
computer system. The bulk of the modeling effort was accom- 
plished on a Hewlett-Packard 9830 desk top calculator with 
associated graphics equipment and the PSAP1 graphics package 
[Ref. 2] implemented on the NPS computer by Lt. Adrian Kibler. 
Material properties were those of hot-pressed Silicon-Nitride 


provided by AIRESEARCH. Post-processing graphics utilized a 


a 


contour plotting routine developed by Gary L. Giles of Langley 
Research Center [Ref. 3] and implemented on the NPS computer 
by Dr. Gilles Cantin of NPS and the author. 

The complete analytical efforts were divided into four 
major categories: 

1) development of a discretized model of the blade and 
attachment design. 

2) development of a pre-processor to calculate consist- 
ent centrifugal loads. 

3) modification of the ADINA code to yield a maximum of 
stress output locations. 

4) execution of ADINA with model developed and elab- 


oration of post-processors to help in the interpretation of 


results. 


II. ANALYTICAL MODEL DEVELOPMENT 


Developing the finite element mesh of the complex geo- 
metric shapes incorporated in the blade was a compromise 
between geometric accuracy, mathematical compatibility and 
economy of effort. Using drawings provided by AIRESEARCH 
[Refs. 4, 5, 6], a mathematical definition of the airfoil 
and the root was developed separately and then mated using 
twenty node bricks arranged compatibly throughout the 


entire assembled mesh. 


A. BLADE NOMENCLATURE 

Figure 1 illustrates the nomenclature used in this paper 
to describe the gas turbine blade model and the reference 
system used to define the geometry. The portion of the com- 
ponent aerodynamically designed to be in the gas flow is 
termed the airfoil. The fillet is the transition section in 
between the airfoil and the attachment root designed to 
minimize stress concentrations. The attachment root (also 
termed the dovetail) forms the base of the airfoil, provides 
radial and tangential placement of the blade, and transmits 
the forces developed by the gas flow to the disk. A right- 
hand cartesian coordinate system was used to describe the 
geometry. The X-axis coincided with the turbine axis-of- 
rotation. The Z-axis coincided with the stacking axis for 
the airfoil profiles. The Y-axis was then mutually perpen- 
dicular to the X- and Z-axis, with the system origin on the 


rotation axis. 
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B. GEOMETRIC DEFINITION 


1. Airfoil 

Reference 4 provided sixty-two perimeter points for 
x-y plane cross-sections at ten z-levels from z=3.1 to z=4.0 
inches. Linear fairing was prescribed for machining between 
cross-sections in order to facilitate manufacture. The 
airfoil tip level of 3.874 inches and z-level 3.385 inches 
(level of beginning of fillet definition) were plotted and 
stored for use in discretizing the airfoil. Because of the 
linear fairing technique, intermediary points were left to 
the ADINA code feature of node generation for definition. 

2. ‘prbet 

The fillet geometry consisted of two radii; 0.3 inch 
and 0.1 inch. The locus of centers of the 0.3 inch part of 
the fillet was at a constant level of 3.385 inches. The arc 
subscribed was tangent to the linear airfoil. The locus of 
centers of the 0.1 inch radius arc was centered such that 
the arc subscribed was tangent to both the 0.3 inch radius 
fillet section and the base which was considered to be a 
3.2 inch radius right circular cylinder about the axis of 
rotation. This double radius design was chosen in order to 
reduce the possibility of stress concentrations by approx- 
imating an elliptic geometry. 

Geometric definition of the fillet was made in 
sixty-two, two-dimensional planes parallel to the Z-axis 


and containing the normal vector to each of the six*y-two 


Milde NPR tea Alan ts ? ‘ ‘ coc inated < 


defined perimeter points at 2=3.385 inches (Figure 2a). The 
following described algorithm was developed to define the 
limiting points illustrated by Figure 2b. 
a. Determination of the External Normal 

A second degree Lagrange interpolating polyno- 
mial was passed through three successive perimeter points at 
the cross-section z=3.385 inches, and the first derivative 
was evaluated for the middle point in order to determine the 
slope of the tangent vector at this point. The negative 
reciprocal of the derivative is the slope of the line normal 
to the analysis point. Taking the y-component as -1 and the 
X-component as the derivative (dy/dx), a normal vector was 
defined. This vector was then normalized to give a unit 
vector. The analysis point was tested to determine the 
desired coefficient signs in order to define the normal 
with the direction external to the body of the airfoil. 
This procedure was carried out for each perimeter point, 
thereby defining sixty-two normals from the airfoil. 

b. Locus of 0.3 Inch Radius Centers 

The coordinate system was translated and rotated 
to yield a working plane which was parallel to the Z-axis 
and contained the unit normal vector for the point being 
analyzed. The curvature center was adjusted in this working 
plane at z=3.385 inches such that the 0.3 inch arc was tan- 
gent to the linear airfoil surface. These contact points, 
defining the start of fillet curvature, were transformed to 


the original coordinate system for future use. 


LS 


ec. Locus of 0.1 Inch Radius Centers 


The 0.1 inch radius curvature centers were 

defined in each of the sixty-two working planes by use of a 
Newton iteration scheme to adjust the center of curvature 
such that the subscribed arc was tangent to both the 0.3 
inch are and the top of the root. These points of tan- 
gencies were determined, transformed to the original coor- 
dinate system and stored. 

d. Definition of Perimeter Point at Arbitrary Z-Level 

Using the above determined points (Figure 2b), 

an algorithm was developed which took as input any desired 
z-level in the fillet region, tested for the correct def- 
inition of geometry for that level and yielded an array of 
sixty-two perimeter surface points for that level. 

Information from the above procedure allowed the 
definition of a sixty-two point perimeter at any z-level of 
an airfoil or fillet cross-section parallel to the x-y plane. 

3. Attachment Root 

The attachment root geometry was defined in a ver- 
tical plane containing the stacking axis and perpendicular 
to the horizontal longitudinal centerline of the body of the 
root which has a broach angle of 23° relative to the axis 
of rotation (Figure 3). Nine different geometries con- 
sisting of straight lines and arcs of circles were defined 
by Ref. 6 in this plane for the perimeter of the root. 
These goemetries were mathematically defined, and the limit- 


ing point coordinates of each segment were determined. This 
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shape was projected into the fore and aft faces of the disk. 


Symmetry conditions allowed simple manipulations of coor~ 
dinates in order to define the shape of the entire attach- 
ment root. An algorithm was developed to define 231 points 


on any desired horizontal surface. 


C. ANALYSIS MESH DISCRETIZATION 

Major considerations in the development of the analysis 
model were to define a mesh which adequately described the 
geometry of the blade and yielded sufficient data points 
for meaningful evaluation, without making the mesh so fine 
as to require an inordinate amount of computer time for 
solution. Twenty node bricks were chosen for use throughout 
the mesh because of their capability of accurately defining 
any second order geometric curve (used almost exclusively 
by the designers of the blade) and their isotropic sensi- 
tivity to the applied loads. 

The airfoil was represented by a mesh of twelve elements 
arranged one deep, three high and four long (Figure 4). The 
top two layers defined the linearly faired section, and the 
bottom row defined the fillet region. Choice of nodal point 
coordinates were made with the assistance of plots of hori- 
zontal cross-sections of points defined by geometry defini- 
tion algorithms. 

Chosen for the attachment root representation were nominal 
three deep, four long layers of elements defining each dif- 


ferent vertical geometric segment (Figure 5). The exact 
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make-up of the root mesh was controlled by the mating of the 


airfoil with the attachment root which caused significant 
distortions because the design definition of the fillet 
base extended beyond the limits of the top surface of the 
attachment root. This mis-match of geometries necessitated 
brute force adjustment of nodal point coordinates in order 
to mate the two regions into a single consistent mesh. 

Once the mating process was accomplished, all the chosen 
nodal points and the element connectivities were card punched 
in the format required by ADINA. The ADINA deck was then 
used as input to the PSAP1 graphics package in order to 
de-bug the mesh. After adjustments were made to yield a 
visually satisfactory mesh, a trial analysis was run which 
brought to light numerous mathematical difficulties caused 
by the distortions resultirg fromthe mating of dissimilar 
geometries. These were corrected, and the final mesh 


(Figure 6) was considered ready for analysis. 
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Figure 1. Illustration of Nomenclature and 
Reference System of Turbine Rotor Blade. 
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Definition of nomenclature for Figure 2 
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outward pointing normal 

level z = 3.385 inches 

center of 0.3 inch fillet radius curvature 

point of tangency between linear airfoil and fillet 
curvature 

center of 0.1 inch fillet radius curvature 

point of tangency between 0.3 inch fillet curvature 
and 0.1 inch fillet curvature 

point of tangency between 0.1 inch fillet curvature 


and attachment root surface 


Figure 2. Continued, Page 2 of 2. 
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(a) Elements 13 through 102 


Attachment Root Mesh, 


Page 1 of 2. 


Figure 5. 
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EASTERLING-FINITE ELEMENT ANALYSIS 
OF A CERAMIC GAS TURBINE BLADE 


Figure 6. Completely Assembled Finite Element 
Analysis Mesh of Proposed Ceramic Gas 
Turbine Blade Design. 


III. PRE-PROCESSOR DEVELOPMENT 


The ADINA finite element code is a flexible, state-of- 
the-art code encompassing linear and non-linear, static 
and dynamic analyses with a choice of six types of elements 
and twenty material models in various combinations with 
element types. In this section, an explanation of the facets 
and capabilities of ADINA used in the blade analysis is 
given along with a description of the preparations that were 


necessary prior to execution of an analysis. 


A. GENERAL DESCRIPTION OF ADINA 
The brittle behavior of the hot-pressed Silicon-Nitride 
lends itself to the use of the linear isotropic elastic 
material model. Twenty node isoparametric bricks were used 
throughout the mesh. ADINA's use of an out-of-core solution 
scheme and compacted storage of matrices lends itself attrac- 
tively for analysis of a large mesh which is required for 
the blade analysis because of the necessity of describing 
many different geometries. 
For the chosen material model and static analysis, the 
equilibrium equation solved by ADINA is: 
KU = R 
where 
K = stiffness matrix 
U = displacement matrix 


R = load matrix 


Pa 


The stiffness matrix is defined as: 


: T 
K = Jie CBdV 


B = strain-displacement transformation matrix 


where 


C = constituative stress-strain matrix 
The elements of B are functions of the natural coordinates 
r, s, t, being derived from the isoparametric representation 
of displacements and the inverse of the Jacobian matrix. 
The integration is carried out in the natural coordinate 


system of reference, and dV is defined as 


dV = det(J) dr ds dt, 


where det(J) is the determinant of the Jacobian matrix. 
Integration is accomplished using Gauss quadrature with the 
resulting matrix elements being stored in compacted form. 

With the exception of gravity loading, concentrated 
nodal forces are the only allowed force inputs incorporated 
in the ADINA code implemented on the NPS computer. The 
rotation induced centrifugal body force needed for this 
analysis was obtained using an external pre-processor to 
be discussed in Section III-C. 

In ADINA, boundary conditions constraining the structure 
model are imposed by eliminating global degrees of freedom 
of selected nodes. The real structure was constrained by a 
Waspaloy disk with a compliant layer interface between the 


ceramic contact surface and metal disk. Interpretation of 
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the behavior of this system of materials varies among inves- 


tigators with consequent variations in modeling the boundary 
conditions. Without concise information on the behavior of 
the actual constraining mechanisms and having to limit the 
scope of this investigation, the tack taken was to eliminate 
all degrees of freedom for all nodes defining the blade 


contact surface. 


B. MODIFICATIONS TO ADINA 

The ADINA code's complex dynamic dimensioning scheme 
inhibits the user from making major adjustments, additions 
or deletions to the implemented code. Two modifications, 
however, were deemed important enough to implement within 
the code. Other problem related manipulations were imple- 
mented by external algorithms which use the ADINA input deck 
but do not affect the code. 

Originally, ADINA came to a STOP statement when a zero 
or negative determinant of the Jacobian matrix was encountered. 
This is a valid procedure since the results of subsequent 
analysis would be incorrect; however, given the inexperience 
of this investigator and the distortions of various elements 
due to the requirement of producing a compatible and complete 
model, numerous elements were suspected of being ill-behaved, 
and the STOP procedure prevented efficient troubleshooting. 
The diagnostic output also inhibited efforts to correct model 
inadequacies because of a lack of identification of the 


offending node and element. These problems were circumvented 
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by deletion of the STOP statement and addition to the diag- 
nostic the location of the offending node and element, 
allowing identification of all improper elements on a single 
trial run. 

ADINA stress output consisted of three normal and three 
shear stresses for sixteen nodes out of a possible twenty- 
seven locations for each element being transmitted solely 
to the line-printer device. The necessity for meaningful 
manipulation of the massive output derived from the 682 
node mesh and the desire for higher density value coverage 
for input into contour plotting post-processing dictated 
modification of the ADINA output capabilities. An addi- 
tional input parameter was defined in the master control 
cards which allowed the user the option of requesting 
output values of all twenty-seven possible output locations 
(Figure 7). A WRITE statement was inserted to output, on 
a user defined file 57, all nodal stresses. Appendix A 
details the changes necessary to the ADINA User's Manual 
[Ref. 1] resulting from these modifications along with other 


changes from a companion investigator's contributions. 


C. CENTRIFUGAL LOADING 

Principal modes of loading considered for the blade were 
thermal, pressure and centrifugal. Given the temperature 
distribution, ADINA has the capacity for determining thermal 


gradient induced loads. External pre-processing was chosen 


| 
| 
| 
| 


for development of pressure and centrifugal loading. The 
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pressure loading pre-processor was developed by Lieutenant 


J. Preisel and reported in Ref. 7. A centrifugal loading 
algorithm was developed by this author and is described 
herein. 

The ADINA version implemented on the Naval Postgraduate 
School computer accepts only concentrated nodal forces. An 
algorithm was developed which would calculate the nodal 
consistent loads equivalent to a rotation induced body force 
and produce an output of punched cards in a format accept- 
able to ADINA. 

1. Mathematical Formulation 

Starting with Newton's Second Law: 


gx L270 


Se oe 


where 
F = external forces 
I = inertia forces 
m = mass 
op = mass density 


a = acceleration 
dv = differential volume 
A vector of discretized consisten nodal forces (V3) is 


desired such that 


<v,>f{u,} = <I>{u} 
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where 


ag = vector 


fu, } = vector 
<I> = vector 
{u} = vector 

<u, V, 


of consistent loads 
of discrete displacements 


of inertia force components I,, Ty. I, 
of displacement function components 


wot 


Using a right-hand cartesian coordinate system and isopara- 


matric element definitions, displacements and coordinates 


are discretized by 


ful} = <hi> tu, 
i i 


{x} 


where 


ai hes = vector 


<h,> {x 


} 


(3) 
} 


of shape functions = f(r,s,t) 


¥ 2 z 
{x} = continuum coordinates <x,y,z> 


{x,} 


discrete nodal coordinates 


The acceleration matrix is defined by 


a=w = Cw x 


where 


In |e 
“ “ 


|m 
u 


Yr) (4) 


angular velocity vector = w,i + Wyd * WK 
: a 


x 


position vector = xi + yj + zk 


c ion = i + 5 + 
acceleration vector ai 20d. alk 


After vector algebra manipulation, one derives 


{a} = [A] {x"} 
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=s< > 
{a} ays ays a, 


{x#}= <x, y; a>? 


Discretizing equations (1) and (5), using equation (3), and 


substituting into equation (2), one has 


“Ni5> {u,} = ofex,> {h;} [A] <h;> dVol {u,} 


dVol = det J dr ds dt 
Simplifying and rearranging yield the following equation 
which was used to program the computer algorithm for develop- 


ment of consistent loads: 


{v.} = oftn CAT <n,» O43 dVol 


Computer Algorithm for Program Centrifugal Load 


The centrifugal load algorithm is generalized for use 
by any ADINA three-dimensional brick element mesh defined 
by eight to twenty nodes per element. A dynamic dimension- 
ing scheme is used in order to conserve core and simplify 


changes dictated by the size of the user's mesh. 


a. Input 

Angular velocity about the three cartesian axes, 
mass density and an ADINA input deck are the required input 
for calculations in Program Centrifugal Load. In addition, 
several housekeeping parameters are required as defined in 
Appendix C. 

The ADINA input deck provides node coordinates 
and mesh connectivity. By using the subroutine RDADIN, 
duplication of the considerable effort to punch and debug 
a large mesh is obviated. This subroutine should prove 
very useful for future users of ADINA who wish to manipulate 
mesh values. 

b. Calculation of Consistent Nodal Loads 

The user can select two to six Gauss point inte- 
gration in the calculations of the consistent loads. Shape 
functions were derived from Ref. 8, which are the same used 
in ADINA. After calculation of loads for each element, the 
contributions to each node are summed in each of the three 
cartesian coordinate directions to yield the desired output 
of concentrated nodal forces. ADINA has the capability of 
summing concentrated nodal forces in a common direction for 
a given node; therefore, loads other than centrifugally 
induced may be input by the user without necessitating mod-~ 
ification to Program Centrifugal Load output. 

ec. Output 

Four options allow the user to restrict the out- 

put to what is desired for the particular problem under 


investigation: 
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a = —aee os a 


(1) ICHK1. This flag allows the option for 
printout of the nodal coordinate data and mesh connectivity. 
Since this data is also output by use of the ADINA code and 
PSAP1, the user would normally not desire this option from 
Program Centrifugal Load. 

(2) ICHK2. This parameter controls the option 
for controlling the output of consistent loads calculated 
for each element. The information may be useful in comparing 
the contributions of various elements. 

(3) ICHK3. The information required by ADINA 
is controlled by this parameter which allows printout of the 
totaled consistent loads for each node for the three coordi- 
nate directions. 

(4) ICHK4¥. This parameter gives the user the 
option of not punching the data cards as may be the choice 
for a first run. 

Other output items include the total force in each 
of the three coordinate directions and the total number of 
load cards punched. The latter value is required by ADINA 
as parameter NLOAD on the Load Control Card described in 


Section IV of Ref. l. 
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Figure 7. 


Stress Output Locations for ADINA 
Linear Static Analysis Results. 
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IV. PROBLEM SOLUTION 


The ADINA problem solution was executed for two, three 
and four Gauss points used for development of the stiffness 
matrix. Static linear analysis was chosen for the solution 
method, and a material model for isotropic behavior was 
used. A Young's modulus of 4.0 x to" psi was used, cor- 
responding to the approximate value for hot-pressed Silicon- 
Nitride at 2500°F. Poissons ratio was chosen to be 0.33. 
Loading was accomplished using concentrated nodal forces 
with values obtained from the centrifugal load pre-processor 
for an angular velocity of 42,000 rpm about the X-axis. 


Table 1 summarizes the execution times for the three 


analyses. 
OPERATION ORDER OF INTEGRATION 

2 3 4 
INPUT PHASE 22.7% 20.83 24.52 
MATRICES 
ASSEMBLAGE (17 oL3 EWS. 63 3305.98 
TRIANGULARI ZATION 
OF STIFFNESS MATRIX 63:7 «36 663.63 105 69:8 
STEP BY STEP 
SOLUTION LO. 8 Peo Nap ee Woe 
TOTAL SOLUTION 
TIME 1549.57 2234.96 4148.07 


TABLE 1. ADINA PROBLEM SOLUTION TIMES (seconds) 
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V. POST-PROCESSOR DEVELOPMENT 


Static linear analysis using the modified version of 
ADINA yields three orthogonal displacements for each 
input node and three normal and three shear stresses at 
each of twenty-seven locations in a twenty node brick. | 
The value of the stress components calculated at a given 
node will, in general, differ from element to element. 
Experience has shown that the straight average of these 
values give the best estimate of the real stress at that 
node [Ref. 8]. For problems involving brittle material, 
the principal stresses (particularly the maximum prin- 
cipal stresses) are the primary values of concern in 
predicting failure of the modeled component. These 
factors necessitated the development of post-processors 
to manipulate the ADINA output into the values required 
for analysis of results. The consequent output of these 
post-processors was still unmanageable from the standpoint 
of visual perception of the relation of the results with 
model geometry. Additional post-processors were, there- 
fore, developed to yield the coordinates and connectivity 
of a twenty-seven node brick mesh and output of punched 


cards in a format applicable to the contour plotting code 


implemented on the Naval Postgraduate School computer. 


A. PROGRAM KCONT 


Prior to data manipulation, the coordinates of all nodes 
of the twenty-seven node brick mesh and the connectivity 
must be known. Program KCONT was developed to accomplish 
this task by using the basic isoparametric finite element 


model equalities: 


n 
x = » a h; Xs 
ea 
n 
: > Ry Yi 
sles 
n 
Zz = = h; Zs 
pak 


where 


X,; Y, Z = global coordinates of any point 
within the element 


Bes Ves 2 = global coordinates of the nodes 
defining the element 


h. = the array of shape factors as 
the function of the natural 
coordinates r, s, t. 

Program KCONT first reads an ADINA input deck. It takes 
the input node coordinates and mesh connectivity data and 
computes the coordinates of the additional nodes. The 
twenty-seven node brick mesh connectivity is assembled and 
is output along with the coordinates of all nodes to the 


line printer and to the user defined storage files in 


accordance with the instructions in Appendix D. 
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B. PROGRAM STRESS 


Program Stress reads the file defined by the ADINA code 
user for storage of stress results and the files storing 
the total connectivity and coordinate information, averages 
all values for each given node and determines the principal 
stresses from the average values. The output consists of 
storage of the results on user defined files and a line 
printer presentation. Refer to Appendix E for instructions 


on using Program Stress. 


C. PROGRAM CONTOUR PLOT DATA 

Data required by the two-dimensional iso-line contour 
plotting code consists of housekeeping information in the 
form of defined NAMELISTS, coordinate data of nodes on the 
analysis plane, four-node evowdamens tonsil connectivity and 
input values of information to be plotted. Program Contour 
Plot Data was designed to define the desired analysis plane(s) 
and yield a punched data deck of coordinates and plotting 
value information. A graphics generated plot of node points 
and global numbers is also generated to assist the user in 
the formulation of the connectivity. The connectivity in- 
formation must then be punched on data cards by the user 
and inserted into the appropriate data deck location. 

Analysis planes can be defined with any of three options. 
The user can simply input the nodes to be used on data cards 
in the format 1615. Most models, however, have surfaces of 


interest which are parallel to one of the three orthogonal 
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planes defined by the axes of the right-handed cartesian 
coordinate systems. Should this be the case, the user may 
define bounding values of a coordinate for testing the mesh 
for the desired nodes defining the analysis plane. If the 
user inputs both bounding values equal, an equality test 

is made in order to define the desired plane. 

The program as written takes the array of plotting 
values from the file containing the principal stresses 
established by Program Stress and uses the maximum prin- 
cipal stress array (variable SIGMX) for plotting contours. 

A user may redefine the read statement (line CTRP2470 of 
Appendix F) in order to input the desired values for a 
particular problem, filling the array SIGMX with the desired 
plotting data. Coordinate values are read from file 58 
established by Program KCONT. 

The resulting punched deck of cards requires the input 
of connectivity and NAMELISTS €0PTION and €&PICT in locations 


designated by Ref. 3. 


D. PROGRAM CONTOUR PLOT 

The contour plotting routine chosen to display stress 
results was developed by Gary L. Giles of the Langley 
Research Center for use with a variety of graphic systems 
including the CALCOMP plotter installed at the Naval Post- 
graduate School. Unfortunately, the NPS graphics software 
package rotates the plotting axis 90° clockwise from the 


original software conception of plot orientation. In order 
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to achieve optimum utilization of the plotting surface, the 


decision was made to modify the contour plot routine to 
utilize the NPS graphics package but yield a plot oriented 
with the horizontal plotting axis along the length of the 

paper roll as originally conceived by the CALCOMP manufacturer. 
In general this requires inputting to the drawing routines 

of NPS the negative of the desired vertical coordinate in 

the calling location for the x-coordinate and the positive 
horizontal coordinate in the calling position of the y- 
coordinate. 

The plotting origin for the analysis plane is chosen to 
be the geometric center. In order to place the plot properly 
on the plotting surface for any set of coordinates, para- 
meters PXORGN and PYORGN were added to NAMELIST &S0PTION with 
default values of 0.0. Inputting the user's coordinates 
for the center of his analysis plane properly centers the 


plot on the plotting surface. 
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VI. RESULTS OF ADINA STRESS ANALYSIS 


Stress results were analyzed using a program developed 
to compare the results of the two, three and four point 
Gauss integration analyses and plotting the iso-maximum 
principal stress contours for various chosen planes of 
analysis. Severe surface tensile stress concentrations 
were observed in the region immediately above the ceramic- 


disk contact region. 


A. COMPARISON OF MAXIMUM PRINCIPAL STRESSES 

Formulation of the stiffness matrix using two-point in- 
tegration resulted in higher stresses caused by a lower 
order in integration, creating a more flexible system which 
concurs with finite element research [Ref. 8]. Some inves- 
tigators have found that the results of reduced order in- 
tegration reflect more accurately the actual component 
stresses because the strict mathematical development of 
the finite element method creates a system which is stiffer 
than the actual component. Insufficient experimental data, 
uncertainty of boundary conditions and the lack of an ade- 
quate convergence study prevent a statement of qualitative 
opinion of the actual stresses involved for the blade 
analysis; however, noting the failure rate of the few 
specimens which have been tested, it is probable that 


tensile stresses in the dovetail are higher than predicted. 


Table 2 illustrates the differences between the three 


analyses by comparing the greatest maximum principal stress 
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encountered, the average of maximum principal stresses and 
the Euclidean norm of the differences between results for 


different integration order analysis. 


ORDER OF INTEGRATION 


2 3 4 
MAXIMUM PRINCIPAL 
STRESS 52862.39 46671.73 4674.73 
AVERAGE MAX 
PRINCIPAL STRESS 6198.79 6100.42 6093.98 
INTEGRATION ORDER ANALYSES COMPARED 
%a3 2-4 3-4 
EUCLIDEAN NORM 122.2901 122.6712 4.7777 


TABLE 2. COMPARISON OF RESULTS OF ANALYSES USING 
TWO, THREE AND FOUR POINT INTEGRATION (psi) 

The maximum value occurred at node 394 which is located 
immediately above the contact region at level 2.9029 (Figure 
8). Results from the three and four order of integration 
analyses showed 1: :ctle difference, indicating that the 
order of integration necessary for the "exact" evaluation 


of the stiffness matrix elements is being approached. 


B. CONTOUR PLOT ANALYSIS 

Thirty-four analysis planes were defined for contour 
plots. Presented herein are eight of these planes which 
show the regions of highest stress concentrations and the 


distribution of stresses in the airfoil. Figure 9 illus- 


trates the relative position of the six plots of analysis 


planes in the attachment root. The level number refers to 


the z-coordinate value of the nodes on that plane. The 
figures presented consist of an element plot, illustrating 
the coverage of input values to the contour plotting code, 
and the iso-stress plots from the second and fourth order 
integration analyses. 

Figure 10 shows severe surface stress concentrations on 
both pressure and suction sides of the dovetail in the 
region immediately above the blade-disk contact surface. 
This result is consistent with other analyses of blades of 
similar design regardless of the boundary conditions used 
and loading scheme, indicating the stress distribution is 
principally a function of geometry. The orders of magnitude 
of the stress concentrations are fifty ksi for the two-point 
integration plot and forty-five ksi for the four-point 
integration plot. 

Figures 1l and 12 present another view of the vertical 


stress distribution in dovetail. Noteworthy information 


from these plots is that the stress concentration is 
generally uniform along the length of the attachment root. 
Figures 13 through 15 are plots of three horizontal 


| surfaces in the region of high stress concentration in the 


: dovetail. 

Another region of concern in the design of a ceramic 

| turbine blade was the fillet area of the airfoil. Figures 
16 and 17 illustrate the stress distributions in the airfoil 


as viewed from the pressure and suction sides. Some stress 
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concentration does occur at the leading and trailing edges 
at the mid-fillet height; however, the values are minor 
compared to magnitudes obtained in the attachment root and 


relative to the strength of the material. Test results to 


date also indicate the adequacy of the fillet design as no 


failures have originated in this region except for impact 
initiated failure caused by flying debris upon failure of 


other blades. 
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Figure 19. Continued, Page 2 of 3. 
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lines for 4-pt Order Integration (ksi) 


Figure 10. Continued, Page 3 of 3. 
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VII. CONCLUSIONS 


The time required for development of pre- and post- 
processors prevented a complete analysis which would have 
the potential for strong conclusive opinions and possible 
suggestions for design improvement. Software and hardware 
tools are, however, now implemented and identified which 
will allow follow-on analysis by NPS thesis students of 
geometric complex machine components. Despite the handicap 
of an inconclusive analysis, certain comments can be made 


concerning the results. 


A. ANALYSIS RESULTS 

Three other finite element analyses of ceramic turbine 
blade designs were reviewed by the author [Refs. 9, 10, 224. 
Despite differing boundary conditions and loading schemes, 
a surface stress concentration persists in the region imme- 
diately above the blade-disk contact surface. This fact 
indicates that further research is required into the possi- 
bility of optimizing the dovetail geometry in order to 
spread the stresses over a greater bulk of material. If 
such an optimized design proves infeasible or non-existent, 
then the conclusion of this analysis and review indicates 
that hot-pressed Silicon-Nitride is unsuitable even for the 
limited goals of the current project. Using the four point 
bending strength of seventy-two ksi [Ref. 9] as a reference 
strength parameter, a simplistic factor of safety of 1.57 


is realized. Considering the analysis did not impose 
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pressure and temperature gradient induced stresses, this 
factor of safety appears to be quite inadequate for pre- 
dicting a decent probability of success. 

Despite somewhat severe distortions in many of the 
analysis model elements, no mathematical singularities 
were encountered in using a reduced integration order. 
Further analysis using finer meshes may show that the 
reduced integration results are more accurate than the 
"exact" integration order. Should this be the case, approx- 
imately fifteen per cent higher values would be predicted 
than the results of "exact" integration of the stiffness 
matrix elements. These higher values would necessarily 
reduce the probability of success calculations normally 


used in design. 


B. OPPORTUNITIES FOR FURTHER RESEARCH 

The computer hardware and software at NPS, together with 
the groundwork laid by this thesis, allows for many avenues 
of follow-on research of which the following is a partial 
Lis es 

1) refinement of the mesh of this analysis and con- 
vergence study of results, 

2) analysis of the blade loaded under pressure and 
temperature gradient induced stresses, 

3) optimization of attachment root design in order to 


alleviate stress concentrations, 
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4) probabilistic failure analysis using ADINA generated 
stresses, 

5) frequency analysis of the blade, 

6) investigation of appropriate boundary conditions, 

7) analysis of the three material system of the blade- 


disk contact region. 
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APPENDIX A 
ADDENDUM TO ADINA USER'S MANUAL, REPORT 
82448=-1, MIT, SEPTEMBER 1975 
(Revised May 1976) 
Two additional output options are available to the ADINA 
user at the Naval Postgraduate School besides what is 
described in the published user's manual. Master Control 


Card one, documented on page II.1, may be altered as 


follows: 


NOTES COLUMNS VARIABLE ENTRY 

(a) 71-75 LTPS 7 Indicator for extended stress 
output and storage of stresses 
on user defined file 57. 
EQ.03; option not desired 


EQ.1; option desired 


(b) 76-80 ITPS8 Indicator for storage of dis- 
placements on user defined 
file 58. 


EQ.0; option not desired 


EQ.1; option desired 


NOTES: 

(a) The use of this option gives the user an output 
of stresses at all twenty-seven locations per element 
described in Section X.3 and illus*rated on page X.37. 


Stresses are stored sequentially on user defined file 57 


in the following order: 


INPUT STEP 


Be element number 

2 integer 1 corresponding to figure X.5 

3. normal stress for node in global X direction 
4, normal stress for node in global Y direction 
5 normal stress for node in global Z direction 
6. shear stress XY 

ie shear stress XZ 

8. shear stress YZ 


Steps 2 through 8 are repeated for the remainder of the 
twenty-seven nodes of the element, followed by a complete 
cycle for each of the remaining elements. 

(b) Use of this option stores the six global displace- 


ments for each node on user defined file 58 in the following 


' 
i 


order: 
INPUT STEP 
Era node number 
Zi X direction displacement 
a. Y direction displacement 
4, Z direction displacement 
5. X-axis rotation (if present) 
6. Y-axis rotation (if present) 
Wes Z-axis rotation (if present) ' 


Steps 1 through 7 are repeated for each of the input 
nodes. Displacements for the additional nodes for which 
stress results are calculated with option ITP57 are not 


calculated. 
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APPENDIX B 


CONVERSION OF CALCOMP PLOTTING ROUTINES 
FOR USE ON A TEKTRONIX 4012 TERMINAL 


The TEKTRONIX 4012 Computer Display Terminal is the 


most powerful interactive graphics system available to the 


general user at the Naval Postgraduate School. The major 


advantages of this unit are its speed and access to the 


IBM 360 computer via the CP/CMS system. A major disad- 


vantage is there is only one terminal designated for the 


general user. During the final stages of preparation of 


this report, PSAP1l and Program Contour Plot were modified 


for use with the TEKTRONIX terminal in order to achieve a 
higher quality of plot in a timely manner than is possible 
with the CALCOMP system. These modifications proved quite 
simple and resulted in a usable plotting system. Presented 
in this appendix are the subroutines added to the graphics 
programs used by the author and a few general comments. 
The user should consult Refs. 12, 13, 14 and 15 for detailed 
information on using the TEKTRONIX unit. 

The user must first determine the plotting origin of 
; his routine and establish the limits of his plotting window 
using subroutine VWINDO. Caution must be used in order to 
ensure that the original program will properly scale the 
plotting values in both coordinate directions. CALCOMP PLOT 
statements which shift the origin must be deleted. Sub- 
routines must then be added to the routine which are titled 


with CALCOMP plotting library routine names which in turn 


call equivalent subroutines in the TEKTRONIX library. The 


following are examples used with PSAP1 and Program Contour 
PLoOes 
1) Subroutine equivalent to CALCOMP PLOTS: 


SUBROUTINE PLOTS 
IBAUD=480 

CALL INITT (IBAUD) 
CALL BINITT 
RETURN 

END 


2) Subroutine equivalent to CALCOMP PLOTE: 


SUBROUTINE PLOTE 
IX=0 

IY=780 

CALL FINLTT CIx,1Y) 
RETURN 

“END 


3) Subroutine equivalent to CALCOMP LINE: 


SUBROUTINE LINE (Y,X,N,N1,N2) 
DIMENSION X(21), Y(21) 

CALL MOVEA (X(1),Y(1)) 

DO 10 I=2,N 

CALL DRAWA (X(I),Y(I)) 
CONTINUE 

RETURN 

END 


4) Subroutine equivalent to CALCOMP NUMBER: 


SUBROUTINE NUMBER (Y,X,H,AL,TH,N) 
DIMENSION IARRAY (4) 

CALL MOVEA (X,Y) 

CALL IFORM (AL,4,IARRAY, 32) 

DO 10 I=1,4 

IF (IARRAY(I).EQ.32) GO TO 10 
CALL ANCHO (IARRAY(I)) 

CONTINUE 

RETURN 

END 


5) Subroutine equivalent to CALCOMP SYMBOL: 


SUBROUTINE SYMBOL (IY,IX,H,BCK,TH,N) 
DIMENSION IARRAY (80), BCK (1) 

CALL CONETA (BCK,N,IARRAY) 

CALL NOTATE (IX,IY,N,IARRAY) 

RETURN 

END 


6) Subroutine equivalent to CALCOMP PLOT: 
SUBROUTINE PLOT (Y,X,I) 
IF (I.EQ.2) CALL DRAWA (X,Y) 
IF (1 .EQ.2) CALL MOVEA (X,Y) 
RETURN 
END 
These subroutines were tailored for use only with PSAP1 
and Program Contour Plot and are much more restrictive than 
the CALCOMP equivalent. For example, subroutine NUMBER 


defined above will only plot integers up to four digits 


whereas subroutine NUMBER used with the CALCOMP will plot 


any integer or floating point number input by the user. 


The subroutines do, however, provide an example of the 
simple modifications required and the TEKTRONIX software is 


available for a much higher degree of graphics sophistication. 
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APPENDIX C 


PROGRAM CENTRIFUGAL LOAD LISTING 
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